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A finite-dimensional Markovian open quantum system will undergo quantum jumps between pure 
states, if we can monitor the bath to which it is coupled with sufficient precision. In general these 
jumps, plus the between-jump evolution, create a trajectory which passes through infinitely many 
different pure states, even for ergodic systems. However, as shown recently by us [Phys. Rev. Lett. 
106, 020406 (2011)], it is possible to construct adaptive monitorings which restrict the system to 
jumping between a finite number of states. That is, it is possible to track the system using a finite 
state machine as the apparatus. In this paper we consider the question of the stability of these 
monitoring schemes. Restricting to cyclic jumps for a qubit, we give a strong analytical argument 
that these schemes are always stable, and supporting analytical and numerical evidence for the 
example of resonance fluorescence. This example also enables us to explore a range of behaviors in 
the evolution of individual trajectories, for several different monitoring schemes. 
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I. INTRODUCTION 



An open system is one which continuously exchanges 
information with its environment [1- 3 . For a Markovian 
(memoryless) system, this amounts to a loss of informa- 
tion from the system into the environment. By moni- 
toring the environment, it is possible to regain this lost 
information and hence to track the system. If the moni- 
toring process is perfect, then one can expect to perfectly 
track the system; that is, to know as much as is possible 
to know about the system, as would be the case it were 
a closed system prepared at will. 

These very general considerations apply equally to 
classical and quantum systems. However, there are some 
very significant differences between the two cases. In the 
classical case, there is only one (best) way to monitor 
the environment. Also, if the classical system has only 
finitely many possible states — this is known as a finite 
state machine [4 — then obviously it is possible to keep 
track of its state using only a finite classical memory — 
another finite state machine of the same dimension. In 
the quantum case, by contrast, there are infinitely many 
inequivalent ways of monitoring the environment that en- 
able the experimenter to deduce what pure state the sys- 
tem is in [H 131 H]. This is because of the entanglement 
between the system and its environment. But in almost 
all cases, an infinite classical memory is required to store 
that pure state, even for a finite-dimensional quantum 
system. 

This last point can perhaps only be understood by 
introducing a little formalism. We consider finite- 
dimensional systems that undergo Markovian open quan- 
tum system dynamics, described by a Lindblad-form 
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master equation [H |3] : 

L 

p = Cp= -i[HeEp - pHl^] ^^cipc\, (1) 

where ^eff = H — i Xlz ^l^i/"^- Here H is Hermitian (it is 
the Hamiltonian) but the jump operators {q : /} are com- 
pletely arbitrary. This decomposition defines the evolu- 
tion of the system conditioned on a perfect monitoring of 
the bath quanta. If a quantum in channel / is observed 
at time t^, the system state jumps from the pre-jump 
state \tlj{t~)) to the post-jump state \tlj{tn)) oc Q|?/;(t~)). 
Then until the next jump its subsequent evolution would 
be generated by the effective (non-Hermitian) Hamilto- 
nian i^eff appearing in Eq. ([l]). 

In general, the post-jump state will depend on the pre- 
jump state |?/^(t~)), and it will not remain stationary un- 
til the next jump, unless it happens to be an eigenstate 
of i^eff- It is thus not at all obvious whether for a gen- 
eral open quantum system it would be possible to keep 
track of its pure state, even in principle, with a finite 
classical memory. On the face of it, it would seem nec- 
essary to store the nature and exact times of each jump 
— requiring a sequence of real numbers {tn : n} each of 
which would require, in principle, an infinite memory to 
store. Alternately one could store the conditioned quan- 
tum state |'0(t)) itself, but this (a finite vector of complex 
numbers) would also require an infinite memory. 

In Ref. [6 we showed that, for an arbitrary D- 
dimensional quantum system obeying a Markovian mas- 
ter equation with a unique stationary (mixed) state, 
one can expect there to exist a monitoring such that 
it is possible to track the system's conditional (pure) 
state with a i^-state machine as apparatus, with some 
(1^ — 1)^ + 1. For a qubit we proved that this is in- 
deed always the case: a two-state apparatus can be found 
that ensures the qubit jumps between just two states. 
This apparatus must implement an adaptive monitoring, 
choosing how to measure the environment depending on 
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its own internal state. This result shows that the infinite 
amount of information required to track a quantum sys- 
tem under a generic monitoring scheme arises from the 
poor choice of scheme and is not intrinsic to the coupling 
between the quantum system and its environment. 

The goal of this paper is to investigate the stability of 
such finite-state monitoring and provide details on how 
such schemes can be identified and constructed. The first 
goal is necessary to establish that the schemes introduced 
in Ref. [6^ are not just mathematical constructions, but 
are physically realizable. For this purpose we restrict 
to qubit evolution, and cyclic-jump schemes. All the 
schemes we study are stable in a mean-square sense, and 
we give a strong analytical arguments plus supporting 
numerical evidence (for two- and three-state machines) 
that such schemes are always stable. 

The paper is organized as follows. In section [IT] we 
review the background for this paper and summarize rel- 
evant results from Ref. 
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we develop our 



[6 . In section 

main results for stability. We discuss these results in the 
context of specific example of resonance fiuorescence in 
section[lVl We conclude with a summary and a statement 
of open problems in this area. 



II. THE PREFERRED ENSEMBLE 

The idea of tracking the system with finite classical 
memory is closely related to the concept of a preferred en- 
semble. To explain the latter notion, we consider an evo- 
lution generated by a Markovian master equation with 
unique steady state defined by Cp^^ = and we assume 
that pss is a mixed state. This mixed state can be de- 
composed in terms of an ensemble of pure states \v^) with 
some positive weights pk such that 



K 



(2) 



k=l 



Note that there are infinitely many such decompositions 
as the states \vl) need not be orthogonal. A valid inter- 
pretation for any such decomposition of a mixed state is 
that in the long time limit a measurement performed on 
the environment will collapse the system into one of the 
possible states \vl) and this will happen with probability 
pk. Such a measurement will, in general, require simulta- 
neously measuring parts of environment that interacted 
with the system at different times in the past. 

A different question arises if we require the measure- 
ments on the environment to be continuous in time, so 
that the system state is continually being collapsed, and 
so (in the long time limit) will be in a stochastically evolv- 
ing pure state. Because the system evolution is Marko- 
vian, this can be done while leaving unchanged the av- 
erage evolution of the system ([T]). Now the natural in- 
terpretation for the ensemble {pkA^k)} ^^^^ ^^^^ 
system will only ever be in one of the states and 
will spend a proportion of time in that state equal to pk. 



In this case, this interpretation is not valid for most en- 
sembles. Decompositions {pkA^k)} ^^^^ realized 
through an experiment via continuous measurement are 
called physically realizable (PR). The fact that some en- 
sembles are not PR is known as the preferred ensemble 
fact [7 . We note here that if we know of a PR ensemble 
with K finite, then if at some time t the system is a pure 
state from the PR ensemble, {v^), and is subject to the 
continuous monitoring that realizes this PR ensemble, 
then in its subsequent conditional evolution the system 
will only occupy states from the PR ensemble and we 
can therefore track such evolution with a classical regis- 
ter with only K states. Such a classical device is known 
as a finite state machine. 

From Ref. [7 we know that the ensemble {pk^ \^k)} 
PR if and only if (iff) there exists rates Kjk > such that 

K 

yk,C\vl){vl\ = J2njk {\v]){v]\ - \vl){vl\) , (3) 

k=l 

where the system jumps between K different states. Typ- 
ically, most ensembles {pk^ \^k)} representing pss are not 
PR, including the K = D ensemble composed from the 
diagonal basis for p^s [3 • 

For a qubit, the conditions in Eq. ([3| can be simplified 
further. Using the Bloch representation, Eq. ([T]) becomes 



(4) 



where A is a 3 x 3- matrix that dictates the evolution of the 
state and 6, a 3- vector, determines the steady state of the 
system, f^s — —A~^h. This equation has a unique steady 
state iff (if and only if) the real parts of all eigenvalues of 
A are negative. We can track this system with a i^-state 
memory iff there exists an ensemble {p/e, r/c} and rates 
Kjk ^0 such that 



Tk • n =1 yk 

K 

Afj + ^ = ^ i^jk{rk - rj) Vj. 



(5) 
(6) 



k=i 



These equations are better than conditions from Eq. (|3| 
for numerical search for PR ensembles because they 
generically reduce the number of different equations, the 
number of different variables as well as the maximal de- 
gree of this system of equations. Thus search for PR 
ensembles reduces to finding a real solution to a system 
of quadratic polynomials with real coefficients. Unfor- 
tunately, this is still a hard problem, which in general 
is known to be NP-complete [15]. Equations ([5])-(|6| are 
used in Sec. IV C to analyze PR ensemble for = 3 for 
a qubit. 

For the moment, we concentrate on a = 2 PR 
ensemble. In Ref. we showed that such ensem- 
ble always exists for a qubit. This can be easily de- 
duced from Eqs. ([5|-(|6|, which for = 2 imply that 
A{fi — r2) = {i<ii2 + /^2iT(^2 — ^i)- This is an eigenvalue 
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equation and we can conclude that 

n=^ss + ^i^ (7) 
r2 =rss - (8) 

where u is the normalized eigenvector of matrix A and 
parameters r]i and r]2 relate to probabilities pi and p2 
for occupying state fi and r2, respectively, and pj = 
Vj/iVi + ^2)- Details expressions for 771 and r]2 can be 
found in Ref. but are not relevant to this paper. Here 
we just need the general structure of the solutions. Be- 
cause the Bloch vectors, fi and r2, must be real, only real 
eigenvectors u of A can contribute to the solution. As a 
3x3-matrix, A has three eigenvalues and, by fundamental 
theorem of algebra, at least one eigenvalue (and conse- 
quently one eigenvector) is real. Therefore, a qubit al- 
ways has a preferred ensemble comprising just two states 
and there can be up to three different solutions. 

Different PR ensembles (including three different PR 
ensembles to track a qubit with one bit [6 ) arise from 
the freedom that experimentalists have (in principle) to 
monitor the system's environment in different ways. This 
freedom exists because of the invariance properties of the 
master equation, Eq. ([T]), which is invariant with respect 
to transformations 

{ci} = +/im| (9) 

. M 

H^H' =H-\Y.^lC^'^-liraC'J) (10) 

m—l 

where jl = (/ii, . . . , /iAf) is an arbitrary complex vec- 
tor and S is an arbitrary semi-unitary matrix — 
Sm=i ^I'm^rni = Realizations of this master equa- 
tion with {c^} as the jump operators and H^^^ = H' — 
^ X]m=i as the effective non-Hermit ian Hamilto- 

nian generate the same average evolution as the original 
master equation, but clearly give rise to different stochas- 
tic evolution. To obtain the most general pure-state un- 
ravelling of the master equation, we require jl and S to 
depend upon the previous record of jumps. That is, we 
require an adaptive monitoring. Of course when we use 
this to achieve jumping between a finite number of states, 
the classical i^-state memory that stores which state the 
system is currently in carries all the information neces- 
sary for determining jl and S. That is, the adaptive un- 
ravelling is specified by k different values for jl and S. 

The physical meaning of these parameters is most eas- 
ily explained in a quantum optics context. The matrix S 
describes a linear interferometer taking the field outputs 
from the system as inputs and interferes them prior to 
detection. This is to be understood in the most general 
sense, including frequency shifters if the system has out- 
puts in different frequency bands. The vector jl describes 
adding (weak) local oscillators to the output fields from 
the interferometer prior to detection by photon counting. 
There have been for many years theoretical proposals for 



adaptively controlling the local oscillator amplitude [8^ or 
phase m E] , and more recently a number of experiments 
have been performed p!QHT2] . one of which (Ref. pT] ) 
used a weak local oscillator, with amplitude comparable 
to that of the system. 

One characteristic that sets apart different solutions 
is the Shannon entropy. Under continuous monitoring, 
the system will occupy states \vl) with probabilities pk- 
The Shannon entropy for an ensemble {pk-, \'^k)^^k\} ^^^^ 
represents pss is h{{pk}) = - E/c Pfe log2 Pfe- This is 
lower bounded by the von Neumann entropy of Pss- 

h{{pk}) > S{pss) = -Tr[psslog2pss], (11) 

with equality iff {pfc, |'^^)('^^|} is the diagonal ensemble. 
We showed in [6 , that some of the K = 2 and K = 3 
ensembles for a qubit can have entropy h much smaller 
than 1. In this case one can track the state of the qubit 
with less than one bit on average, meaning that the state 
of N identical qubits subject to the same independent 
monitoring can be tracked with Nh <^ N bits. This is 
one reason why studying different PR ensembles for the 
same system is of great interest. In the next section we 
explore the stability of such adaptive monitoring. 



III. JUMP DYNAMICS AND STABILITY 

As explained above, the existence of a PR ensemble 
{pkA^k){^k\} ensures that if we start the system in state 
\vl) and subject it to the adaptive monitoring determined 
by the parameters tZjk from Eq. ([3|, then the system will 
always jump between states \vl) with k = 1, . . . , and 
throughout its evolution will never leave the PR ensem- 
ble. But what happens if the system is not initialized 
perfectly at the start of monitoring procedure? 

We address this question in this section. We consider 
a general qubit system, whose evolution is governed by 
the master equation, Eq. ([T]), with one jump operator, 
c. We assume that this system is subject to adaptive 
monitoring that allows the system to jump between K 
states in a cyclic manner. Since the resulting evolution 
is stochastic, we cannot say what will happen to each 
specific trajectory, but we can determine what happens 
to many different realizations on average. We show that 
on average any initial state subject to the adaptive mon- 
itoring will eventually converge to the jumping between 
states \vl) with k = from PR ensemble. We 

first develop our results assuming that adaptive monitor- 
ing leaves the qubit jumping between K = 2 states and 
then generalize to iC- state cyclic jumps. 



A. Two-state jumping 

For two- state jumping scenario, the system jumps be- 
tween pure states \vl) and I'l'l). In principle, we can 
compute these states using Eqs. ([7|)-([8|. This approach 
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tends to yield useful numerical answers, but extracting 
simple analytical expression is not easy. An alterna- 
tive approach for identifying jumping states (for the spe- 
cific example of resonance fiuorescence) was undertaken 
in Refs. [13, 14 . This approach uses the properties of 
transformations, Eqs. ([9|-([l0|, that leave master equa- 
tion, eq. ([T]), invariant, to explicitly construct the adap- 
tive monitoring schemes that generates a PR ensemble 
with the desired number of elements. This is the ap- 
proach we use in this section. 

The system we study in this paper is simple. It has 
only one jump operator and, therefore, the only degree 
of freedom for generating Eqs. ([9|-(10) is the strength 
of the local oscillator, a complex scalar fi. For a given /i, 
the effective Hamiltonian for the system is [3] 




l^2i-l> 




FIG. 1: A schematic of the conditioned state evolution under 
two-step adaptive monitoring. 



H{lJi) = H — -c^c — iii*c - 



Inl- 



and the jump operator is 
c'(m) 



(12) 



(13) 



Note that although for convenience we are not using 
the "eff" subscript anymore, ^(/i) is the non-Hermitian 
Hamiltonian f^eff introduced earlier. 

Under two-step adaptive monitoring, the signal from 
the system (a qubit) is mixed with the local oscillator 
with strength fii prior to the photon detection. Before 
the photon is detected the system undergoes smooth evo- 
lution governed by Hi = H^jii). Then the system expe- 
riences a jump governed by Si = c'^iii) when the photon 
is detected. At this point an experimentalist switches 
the strength of the local oscillator from jii to /i2 and the 
consequent smooth evolution is generated by the effective 
Hamiltonian H2 = ^(^2) and the next jump is caused 
by 52 = c'{ii2)- As soon as this jump is detected, the ex- 
perimentalist switches the strength of the local oscillator 
back to jii and the evolution cycle repeats itself. 

We can solve for the evolution of the system as fol- 
lowing. Let the initial state of the system be We 
assume that the nth jump happens at time tn after which 
the system will be in state 1?/^^) and the waiting times be- 
tween jumps are given by Tn = tn — tn-i- The evolution 
is stochastic and for each realization of such evolution 
jumps will happen at different times. The probability 
density function, p(r), for waiting times r = (ri, . . . , r^) 
is given by 



p{t) 



(14) 



Here the unnormalized state \ipn) is defined recursively 

by 



^"IV^n-i) with j 



1 for n odd 

2 for n even 



(15) 



This construction implies that |V^n-i) implicitly depends 
on Ti, . . . , Tn-i and, consequently, state \ipn) depends on 



all prior waiting times r. The diagram for this evolution 
is depicted in Fig. [l] At an arbitrary time t = tn-i + 5, 
with < s < Tn, for some n, the state of the system is 
given by 



\m) 



-iH\s 



iV^n-l) 



-iH\s 



Hn-l)\\ 



with j 



for n odd 
for n even 



(16) 

We can force the system evolution to be restricted to 
just two states, \vl) and I'^f) if we require that the state 
\vl) is an eigenstate of Hi^ the state \v2) is an eigenstate 
of i^2 5 the jump operator si maps \vl) to \v2) and 52 
performs the reverse action. A consequence of these as- 
sumptions is that after two jumps the state \vl) should 
return back to itself, i.e. 



SlS2\Vi) 



[/il/i2 + (M1 



l) ^ K) (17) 



Without loss of generality, we can assume that the jump 
operator c is traceless and this means that the square 
of this operator is proportional to the identity, oc /. 
Thus, Eq. (17) can happen only in two ways: either (i) 
\vl) is an eigenstate of c; or (ii) /ii + /i2 = 0. Case (i) 
holds iff \vl) is an eigenstate of si, which would prevent 
any jumping from happening. Since case (i) is ruled out, 
case (ii) must hold. Thus, /ii = — /i2 and the full cycle 
jump operator S1S2 is proportional to the identity, which 
is a very useful fact for the discussion later on. 

We also comment here how Eqs. ( 14 )-( 16 ) can be used 



for simulation of a piecewise deterministic process and, in 
particular, for two-state jumping. We start the system in 
the normalized state iV^o)- We know that after nth jump 
the normalized state of the system is given by Eq. ( [161 ) 
with 5 = 0. We determine a random waiting time 
according to the cumulative waiting time distribution [2] 



F{Tn) = 1 



^"l^(T„-l))||' 



(18) 



with j = 1 if n is odd and j = 2 if n is even. This 
is done as follows. We draw a random number 77 from 
the uniform distribution over the interval [0, 1] and solve 
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T] = 1 — F{rn) for r^. Between the jumps the state of the 
system is given by Eq. (16). The evolution depicted in 
Figs. [3] and |4] is generated using this procedure. 

The cumulative waiting time distribution, F{rn) can 
be, of course, related back to the probability density p{r). 
Using the induction step 



Pi^n) = P{rn\rn-l) P{rn-l) 



together with the fact that 



d 

dTr, 



F{Tn)p{Tn-l) (19) 



d 

dTr, 



nrn) 



d 



dTn 



lllV^n-l)|P' 



(20) 



we can deduce Eq. (14) 



B. Mean-square stability 

We now know how to describe the state of the system 
subject to two-state monitoring and we know the proba- 
bility density for jumps occurring at times ti, . . . , with 
waiting times r = (ri, . . . , r^). The resulting evolution is 
a stochastic process, as we don't know when the jumps 
occur. Thus it is sensible to look at what happens to the 
system on average, over many different realizations, by 
considering the mean-square stability. 

The mean-square fidelity between the desired system 
state and the actual system state immediately following 
the nth jump, starting from initial state |?/^o) is 



j Crfp{T)\{v'Mn)\ 
rfV|||^„)|| 



2l(«IIV'n)|2 



(21) 



with j = 1 if n is odd and j = 2 if n is even. Thus this 
quantity is easy to compute. At the same time, it is easy 
to see that if 



lim {\m^n)? 



1 with j 



if n is even 
if n is odd 



(22) 



the system will converge on average to a perfect jumping 
scenario. We call this mean-square stability. 

We now investigate under what conditions the system 
has mean-square stability. To do this, we need a way 
to evaluate the dependence of 1?/;^) on the jump times, 
r. We let and \v^j) be the eigenstates of iHj with 
respective eigenvalues A| and A^, where states \v^) are 
part of the PR ensemble^ and states \v^) are ot/ier states. 



not part of the ensemble. It is easiest to compute the 
time dependence of lipn) if we decompose this state in 
terms of eigenstates of Hi for even n (i.e., \vl) and 1'^^^)) 
and in terms of eigenstates of H2 for odd n (i.e., \v2) and 

Thus we let the state of the system after 21 and 2/ + 1 
jumps be respectively 

\^2i) = a2i\vl) ^ P21K) (23) 
\i^2i+i) = a2i+i\vl) ^ p2i+iK)^ (24) 

The dependence on the jump times is hidden in the co- 
efficients a2u <^2^+i7 P2h ci^nd which we determine 
by establishing the relationship between them and their 
dependence on and /^q. To do this, we first define the 
action of jump operators si and 52 is on the basis states 
{\vl), K)} and {|«|), 



h\vl) = 
hK) = 

S2|f|) = 



QliK) 

Qll\vl)+Ql2\V2) 

Qlilvt) 

\v'i) + Ql2K)- 



?21l 



(25) 



The scalar elements Qj^ and can be determined once 

all the parameters of the system, c, and /ii, are spec- 
ified. Eqs. (25) simply state that jumping operators (<si 
and 52) map jumping states {vf and 'Of) to each other and 
the mapping for the other basis states is quite general. 
Using Eqs. (23)-(25) and propagation between states 



given in Eq. (15), we can derive 



P21+1 
hi 



hi-ie' 



Now we can show by induction that 

I 

f^2l=MQl2Ql2yY[^~^^''''~'^~ 



A2r2 7 



(26) 
(27) 



(28) 



P21+1 = /3oQL(Q22Qi2)' n e-^^°^^^e-^?-^^+^ (29) 
j=o 

This information is sufficient to derive the mean-square 
fidehty {\{v]\i>„)\^), using the following useful identity 



IKI^n)|' = 



\Pn\Hi-\o,n (30) 



where j = 1 if n is even and j = 2 if n is odd and 
Oj = {vj\vj) is the overlap function. Now calculating 
the mean-square fidelity becomes easy. According to 
Eq. (14), IllV^n)!!^ is the probability density for jumps 
occurring with waiting times r = (ti, . . . ,rn) and thus 
when integrated with respect to d'^r yields one. Thus 



IKIV'n)^ 



1 - (1 - \Oj 



(31) 
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with j = 1 if n is even and j = 2 if n is odd. Using 
Eq. (I28|), we learn 



Thus the mean-square fidehty is 



|/3orC'. (32) 



\{vt\i,2i)f) = l-\Pof{l-\0^\')C' 



(33) 



and 



Kt'IIV'2/+i)p 



IQ22/^o|'(l-|02nc'. (34) 



Thus coefficient C determines the convergence rate to the 
perfect two-state jumping. 

The absolute value of fidelity is a number between 
and 1 and its average with respect to all possible waiting 
times r is still between and 1. Thus Eq. (33) implies 
that C < 1. Thus unless C = 1, Eq. (|34| implies that 
the mean-square fidelity always converges to 1 and the 
system has mean-square stability. 

We now determine under what conditions parameter C 
is strictly less than one. To answer this question, we need 
one more calculation, which is performed in Appendix A 
and establishes that IQiiQfi P/4ReAf ReAg = 1 and that 
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= Ql^Qh- Thus 



C = 



ReAfReAl 
ReAJReAg ' 



(35) 



In other words, the system will have mean-square stabil- 
ity iff the geometric mean of the real part of the eigen- 
values of the states in the ensemble is smaller than the 
geometric mean of the real part of the eigenvalues of the 
states not in the ensemble. This makes sense, as the 
real parts of these eigenvalues gives the rate at which the 
amplitudes of their respective states decay during the 
between-jump stages of the evolution, as we will explore 



in Sec. HID Thus, the smaller they are, the more stable 
the respective states. As we will see in Sec. |IVB[ the 
condition for mean-square stability can be satisfied even 
though every second stage is unstable. 



C. i^-state jumping 

We now explain how our results can be generalized to 
cyclic i^^- state jumping, as first studied in Ref. pi for 
K > 2. We simply let the index j for expressions in 
sections III A and IIIB range over {!,..., K} instead of 
just taking values 1 and 2. All considerations still hold 
and we can derive Eq. (IsTl) as before, but the index j 



in this expression is now determined according to j = 
{n mod K) + 1. We also need to reevaluate J d'^r\Pn\'^. 
The unnormalized state for the jth step in cycle after n 
repeats is 



with = nK + j and j + 1 stands for (j + 1) mod K. 



Equation (25) becomes 



Sj\Vj) 



Qiil^I+i) 



^2irj + l 



(37) 



and expression for f3nK becomes 



n-l 



PnK = MQI2 ■ ■ ■ Q22r n e-^°"^^+^ • • • e-^-^<^+i)- 
i=o 

(38) 

and thus 



/- 



Q22 ' ' ' Q22 

2^ReAf • • • ReA^ 



(39) 



Thus parameter C that determines stability is redefined 
to be 



C 



Q22 ' ' ' Q22 
2^ReA? • • • ReA^ ' 



(40) 



Just as before we can show that C can be re-expressed 
as 



C 



ReAf • • • ReA^ 
ReAf • • • ReA^ ' 



(41) 



so that C <1 whenever ReAJ • • • ReA^ > ReAf • • • ReA|^. 
The argument proceeds in the same way as for 2-state 
jumping, and the details appear in Appendix A. 

We conjecture that the fact that \vj) form a PR en- 
semble for pss, the unique steady state, ensures that 
^^^k ^ ReA^. Even more generally, we con- 

jecture that every finite PR ensemble is stable. 



D. Stability of an individual trajectory 

It is important to note that mean square stability does 
not imply that in an individual trajectory there will be 
monotonic convergence of the system towards the desired 
states. The fidelity could decrease between the jumps, 
and/or could decrease due to a jump. We begin with the 
first issue. 

We say that the system is piecewise stable if \{vj\tlj{t))\'^ 
is a monotonically increasing function of t during all 
stages tn-i < t < tn^ where j = 1 if n is odd and 
j = 2 if n is even. Otherwise, if the fidelity with the 
desired state decreases between jumps, this constitutes 
an unstable stage in the evolution. As we will show, this 
happens only if n even or n odd, not both. That is, stable 
and unstable stages will alternate. 

The stability of the evolution between jumps is deter- 
mined by the relationship between eigenvalues A| and A^. 
Let us assume without loss of generality that the state of 
the system right after the jump is 



(36) 



(42) 
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Then the unnormahzed state of the system at time r after 
the jump, but before the next one, is 



Again, we ihustrate this for the resonance fluorescence 



example in Sec. IV B Note however that for mean-square 
stable evolution, in the long time limit, converges to 
IVi(r)) = aoe"^^'^|v|) + /3oe"^j"^|vJ). (43) \v'i)^ and so Fi and F2 also converge to one. Thus we 



We now want to know that happens to the fldelity be- 
tween |'0(r)) (normalized state of the system between 
jumps) and the ideal jumping state \v^): 



expect that as the system converges towards the ideal 
ensemble states, observing jumps that decrease fldelity 
becomes less and less likely. 



\{vmr))? = 



m\^{r))? 



Using Eq. (30), we can compute this fldelity to be 



\{vmr))? 



\0i?) 



\aoe 



-^^'>lvp + /3o\v°)\^' 



(44) 



(45) 



Thus we can see that if ReA^ < ReA^ then in the long 
time limit fldelity approaches one. That is, between 
jumps the system will converge to the desired ensem- 
ble state with exponential rate and the fldelity is mono- 
tonically increasing to the maximum value of one. On 
the other hand, if ReA| > ReA^, then the system state 
converges towards the non-ensemble eigenstate, and the 
fldelity converges towards | O j p , which is a quantity less 
than one. Thus the fldelity may decrease during such a 
stage, and we call this an unstable stage in the evolution. 
This is illustrated for the case of resonance fluorescence 
in Sec. ITVBl 

We now address the issue of whether the jump itself 
increases or decrease fldelity. We again introduce some 
notation. Let = ai\vj) + Pi\Vj) be the state of the 

system right before the jump and \<p2) = <^2|'^p + /^2|'^fe) 
be the unnormahzed state of the system right after the 
jump, where k = (j + 1) mod K. We also know that 
14^2) = We want to compare the fldelity right 

before and immediately after the jump and we let 



Fi 



\{v]\4>i)\' 



and F2 



\{VI\4>2)\' 



(46) 



We evaluate \{Vj\(l)i)f and |('^^|02)| using Eq. (30). 
Then fldelity will decrease after the jump, Fi < F2, m 



2)\\' ^ 1^2? I- W 



< 



= B. 



(47) 



IV. RESONANCE FLUORESCENCE 

In this section we apply our results for stability in the 
mean-square sense, and for individual trajectories, to a 
speciflc physical example. In particular, we show that 
there are PR ensembles which are mean-square stable, 
but whose cycles are composed from stable and unstable 
stages. Others are comprised of only stable stages, and 
so the fldelity is piecewise monotonically increasing. It is 
only piecewise because this example also illustrates that 
jumps can decrease the fldelity. 

The example we consider is resonance fluorescence. We 
consider a qubit (for example, a two- level atom) with 
basis states |0) and with a transition frequency uq. 
We assume that qubit is coupled to the continuum of 
electromagnetic radiation and, therefore, decays to |0) 
at rate 7. At the same time, it is driven by a classical 
fleld oscillating at frequency ujq. The strength of the 
driving is quantifled by the Rabi frequency In the 
interaction frame [3^ with respect to the atomic transition 
frequency ujq , the evolution of the oubit is given by the 
master equation of the form of Eq. (Ilj) with H = {Q./2)ax 
and one jump operator c = y^a, where & = |0)(1| and 

o'x = o- -\- . Then matrix A and b in Bloch vector 



equation, Eq. (|4|, are 




and b 



(49) 



The steady state fss = (0, 2jn, -7^)^/ (7^ + 21^^), is a 
mixed state for ^7 7^ 0. The unnormahzed eigenvectors of 
A are ui = (1, 0)^ and u± = (0, 7 ± ^7^ - 161^2, 41])^. 
For e = Q/j and \e\ < 1/4 all three eigenvectors of A are 
real, while for |£:| > 1/4 only ui is real. Note that e here 
is different from e = (1^/7) ^ in Ref. |6]. 



We can easily determine tight bounds for 



VI 1 1^1) IP by observing that 



An 



< 



\s]sj\(pi) < An 



(48) 

where Amin and Amax are the smallest and the largest 
eigenvalues of the operator s^jSj. Thus the lower bound 
of Amin, or upper bound of Amax, is attained if |^i) is 
the corresponding eigenstate of s^jSj. If Amin < 5, it 
is possible to observe jumps that decrease the fldelity. 



A. Two-state jumping 

We now analyze stability properties for two-state 
jumping in this system and we do this with the method 
described in Sec. IIIB and in Refs. |l3l[T4], i.e., we know 



that PR ensemble is constructed from the eigenvectors of 
H{ii)^ or equivalently, of iH{fi)^ where H{/j^) is the effec- 
tive (non-Hermitian) Hamiltonian given by Eq. (12). 

For resonance fluorescence, the operator —iH^ji) has 
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(a) 



(b) 





(c) 



(d) 



FIG. 2: (Color online.) Solid arrows show Bloch vectors for 
2-state jumping. The volume of the sphere at the tip of each 
arrow represents the probability that the qubit occupies the 
corresponding pure state. The dashed arrow is fss- Solutions 
in (a) and (b) arise from ui with e — 1 and e — 0.23, re- 
spectively. This is the solution that exists for all e. Solutions 
depicted in (c) and (d) exist only for e < 0.25 and shown here 
for e = 0.23. The solution in (c) is generated by i2+ and that 
in (d) by U-. 



eigenvectors given by 



|^±(/i))=£|l)+ ± 



1 



|0). (50) 



The corresponding eigenvalues are 



a±(m) 



± 



(51) 



Equivalently, the probability of occupying each state is 
1/2 and Shannon entropy is 1. Thus one bit is sufficient 
to track the state of the system. As shown in Fig. |2|a)- 
(b), the solution is symmetric with respect to x = plane 
and as e increases jumping states move away from each 
other and away from the bottom of the Bloch sphere. 
For jii = Vj^ and jii = V- the jumping states are gen- 
erated by Uj^ and U-^ respectively. Fig. [2|c)-(d) shows 
that such solutions lie in the x = plane. For /ii = z^-, 
jumping states are spread as far apart as possible; the 
system spends most of the time in one state, which is 
nearly aligned with the direction of the steady state, rgg. 
For /ii = z^+, jumping states cluster together around the 
direction of the steady state, r^^\ the probability of oc- 
cupying each state is more similar and becomes equal as 
e approaches 0. Since the probability of occupying each 
state is no longer equal, the Shannon entropy is less than 
one. This means that one could store the state of the 
qubit in less than one bit on average. That is, one could 
keep track of the state of a collection of N identically 
monitored qubit s using only Nh bits. Actually, the en- 
tropy of the solution due to V- is very low. That is, it 
is very close to the von Neumann entropy of the steady 
state, which is a lower bound for the Shannon entropy of 
any ensemble representing for the steady state. For this 
solution, the Shannon entropy in a Taylor expansion in 
terms of e is 

h{u_) = (log2 10-4 log2£:)£'^-48(log2£) £^+0(6:^) (53) 
compared to 

5(pss) = (Iog2l0-41og2£)e'+6(log2£)£'+0(£^). (54) 
B. Stability of two-state jumping 



where here, and in the remainder of the paper, we have 
set 7 = 1 for simplicity. Whether |'u+(/i)) or \v-{jj.)) is 
part of the PR ensemble depends on a particular value of 
/ii. We know already from Sec. |III A| that fi2 = — Mi, so 
without loss of generality we can choose jui to have a non- 
negative real part. Once we impose additional conditions 
from Sec. |IVC| associated with jump operators Si and 
S2, we learn that /ii only assumes values from the set 
{1/2, where 



. VTTvT^16? 



(52) 



The last two values in the set contribute to solutions only 
for s < 1/4. That is, one only obtains a PR ensemble 
for values of jui that are either real or purely imaginary. 
Detailed derivation of these conditions can be found in 
Refs. [131 [14], ^iid the Bloch vectors for all three ensem- 
bles are shown in Fig. [2] 

For /ii = 1/2, the jumping states (in Eqs. ([7|)-([8|) are 
spanned from eigenvector ui and, as explained in Ref. [6 , 
the system spends equal amount of time in each state. 



For /ii = 1/2, the PR states are \v-^{fii)) and 
\v-^{—jiii)) and the full set of eigenvalues from Eq. ( [sT] ) 
can be simplified to 



A, 



4 



- ± — and A_ 
8 2 



4 



8^2- 



Thus ReXf = ReA| = 1/8, while ReAf = ReAg = 5/8. 
Therefore C = 1/25, and the system is mean-square sta- 
ble. Moreover, it is piecewise deterministically stable 
(that is, the fidelity increases monotonically except per- 
haps at jumps.) 

For /ii = the situation is quite different. The 
PR ensemble is comprised from |'u+(/ii)) and \v-{—fii)) . 
It is still mean-square stable, but the short-lived state 
I'^l) = \v-{—fii)) has an eigenvalue such that ReA^ > 
ReA2, where \v2) = |'y+(— /ii)). That is, every second 
stage is unstable. We illustrate this situation in figures [3] 
and [U Here we assume that evolution starts with the 
unstable step; that is, the initial strength of the local 
oscillator is — /ii. We also let £ = 0.1 and the initial 
state is \ipo) = -1.09|v_ (-/ii)) + 0.5|v+(-/ii)). Figureji] 
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FIG. 3: (Color online.) Alternating stages of stable and un- 
stable two-jump evolution, exaggerated by an atypical trajec- 
tory. The blue (darker) line represent the fidelity of the state 
with the ideal PR ensemble and green (lighter) line shows 
when jumps happen (photon count increment dN — 1). This 
is a highly improbable evolution, as the systems spends an 
uncharacteristically long time in an unstable step and unchar- 
acteristically short time in a stable step. Fidelity is decreasing 
during the unstable step. At the second jump there is also a 
decrease in fidelity due to the jump itself. Time t is measured 
in units of 7"^. 



shows highly untypical evolution, where the short-lived 
state |'U+(— /ii)) exists for an improbably long time and 
the long-lived state \v-{—iii)) lasts an improbably short 
time. One can see here that for stages governed by — /ii 
the fidelity decays as the system moves away from the 
desired PR ensemble state. Figure [4] shows typical evolu- 
tion under the same conditions. Here the instability can 
be seen only during the first step, which is short (at time 
Ti ^ l.l7~^). The next jump does not happen for long 
time ( T2 ^ 792.57"^) and it is very short (it lasts for 
r3 ~ 0.87"^). Thus the unstable stages contribute very 
little to the overall evolution of the fidelity. As a result 
it does not disturb mean-square stability for the system. 

For 111 = Vj^ the situation is even more curious. For 
e < Sq ^ 0.243, the PR ensemble is comprised from 
|'U+(/ii)) and |'u+(— /ii)) and for £ > sq, the PR ensemble 
is generated from |'u+(/ii)) and \v-{—iii)) . The system 
is always mean-square stable, and is piecewise determin- 
istically stable for e < eq. However for e > Sq the short- 
lived stage associated with \v-{—iii)) becomes unstable. 
At 6 = A+(— /ii) = A_(— /ii). In Fig. [sj we show 
the entropy for three different solutions with two-state 
jumping and the type of stability they exhibit. 

Recall from Eq. (32) that the deviation of the fidelity 
from unity decays as C\ after 21 jumps. The critical 
constant for all three values of jii can be written as 



C 



iMil 



4ReA?ReA^ 



(56) 



FIG. 4: (Color online.) As in Fig. [s] but showing a typical 
trajectory. The fidelity decreases visibly only during the first 
stage, as stable stages last much longer than unstable stages. 




0.00 0.05 0.10 0.15 0.20 0.25 
s 

FIG. 5: (Color online.) Ensemble Shannon entropy h for the 
three different two- state jumping solutions. A solid line indi- 
cate mean-square stability as well as piecewise deterministic 
stability. A dashed line indicates that the solution has only 
mean-square stability, with one of its stages being unstable. 
The uppermost (yellow) line shows entropy for PR ensem- 
ble with 111 = 1/2. The middle line (blue) describe entropy 
for solutions with /xi = z/+ and the lowest line (red) is for 



clear (although not perfect) correlation between the low 
entropy solutions and the most quickly converging (in 
the mean-square sense) solutions. Note however that the 
most rapidly converging solution is typically not piece- 
wise deterministically stable. 



We note here that from Eq. ( 33 ) , one can see that con- 



We show the convergence coefficient C as a function of 
e for all three possible values of jii in Fig in There IS a 



vergence coefficient C characterizes only certain aspect of 
mean-square stability. In particular, it determines how 
many stages are needed for convergence to the PR en- 
semble. Thus the smaller C in Fig. |6] corresponds to 
fewer stages needed to achieve certain level of conver- 
gence in the mean-square sense. This plot, however, does 
not tell us about the time needed to achieve convergence. 
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FIG. 6: (Color online.) The convergence coefficient C as 
function of e for the three different values for /ii . The yellow 
(straight) line corresponds to /xi = 1/2, the red (lowest) line 
to 111 — V- and the blue (intermediate) line to \i\ — z/+. 
Note that the lowest entropy solutions, which is not piece wise 
deterministically stable, is the most stable (smallest C) in the 
mean-square sense. 



FIG. 7: (Color online.) The asymptotic convergence rate R 
as function of e for the three different values for /ii . The yel- 
low (straight) line corresponds to \i\ — 1/2, the red (lowest) 
line to \i\ — V- and the blue (intermediate) line to \i\ — v^. 
Note that, a high value of R indicates faster convergence (per 
unit time), so this measure of convergence reverses the order 
of the three schemes relative to that in Fig. [6] 



For example, although the solution with /ii = v_ has a 
much smaller value for C than \i\ = 1/2 solution, and so 
requires more jumps to converge, it could have shorter 
stages between every jump so that it approaches fidelity 
in a shorter time than \i\ = v_ solution. To show that 
this is a plausible scenario we consider an asymptotic rate 
of convergence to the PR ensemble, which we define as 
R = — ln(C)/(tas), where (tas) is the expected value for 
the duration of the complete cycle of evolution for the 
corresponding PR ensemble. Because asymptotically the 
actual ensemble converges to the PR ensemble, and using 
the law of large numbers, the number of cycles under- 
gone in the limit t ^ oo converges to t/(tas). Thus from 
Eq. (33), in the long-time limit, the fidelity F(t) between 
the actual conditioned state and the record-determined 
PR state at time t behaves as 

{Fit)) ~ 1 - |/3o|'(l - |Oinexp(-i?t) (57) 

For two-state jumping this asymptotic convergence 
rate is 



i?= -ln(C)/ [(2ReA^)-i + (2ReA|)-i] 



(58) 



and we plot it in Fig. [7| as a function of e for all 
three values of /ii. Here the solution with /ii = 1/2 
has the highest asymptotic convergence rate, with R = 
ln(5)/4, while the solution /ii = has the lowest, 
with ^ as s ^ 0. Specifically, for fii = al- 
though log(C) = 0(log(£2)) > 1, ReAf = 0{e^) so that 
R = 0{e^\og{£^)) 0. This is a complete reversal from 
the results reported in Fig. [6j based on C. 

The vast difference in rates in the limit £ ^ can 
be understood from the nature of the ensembles as il- 
lustrated in Fig. [2] In all cases the the average pss dif- 
fers from |0)(0| only at 0(e), and differs from a pure 
state only at 0{e^), as reflected in Eq. (54). In the case 



/i = 1/2, the ensemble comprises a pair of Bloch vec- 
tors located near f^s in a symmetric fashion in the plane 
perpendicular to x = plane. For jui = Vj^ the Bloch 
vectors lie in the x = plane around rgg in a nearly 
symmetric fashion. In both cases, the system spends the 
same time in each state and thus the expected time for 
the duration of each cycle is independent of e, so the 
rate of convergence is non-zero. For fii = z^-, one of the 
jumping states is nearly aligned with pss [i-e. it has the 
form |0) +0(£)|1)] and the other state approaches |1) as 
6^0. For small the system spends almost all of the 
time in the first state. It jumps to the excited state with 
a rate 0{£^)^ and jumps back with a finite rate. The for- 
mer process is the rate-limiting step, so the system has a 
cycle whose duration tends to infinity as £ ^ 0, giving a 
convergence rate near 0. 

However, we need to be careful in interpreting the re- 
sults from Fig. [7| as it only shows the asymptotic rate. 
Most of the convergence happens in the initial evolution 
— during the first few stages, while the state of the sys- 
tem is drastically different from the PR ensemble. The 
asymptotic rate becomes relevant only when the state of 
the system closely resembles the PR ensemble and, in a 
way, the convergence has already occurred. We, there- 
fore, also consider another measure for the rate of con- 
vergence: Ri = — ln(C)/(Ti). Here (Ti) is the expected 
time for the duration of the first cycle. For two-state 
jumping, this quantity was investigated numerically. 

Our numerics reveal that the initial convergence rate 
Ri shows strong dependence on the initial condition. For 
some initial conditions, Ri for the /ii = solution is 
above that for the jui = 1/2 solution for all e < 1/4. 
Because the and solutions coincide for e = 1/4, 
this implies that, for these initial conditions, and for s 
sufficiently close to 1/4, the /ii = 1/2 solution has the 
slowest initial rate of convergence for some range of e 
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FIG. 8: (Color online.) Quantities relevant to a fidelity de- 
crease upon jumping in the two- state jumping solution with 
/ii = z/_. The blue (lowest) line shows the achievable lower 
bound for the relative norm 1 1 102) 1 1^/| | |0i) 1 The cyan (up- 
permost) line shows the bound B the jump generated by si 
and purple (middle) line is the same bound for jumps gener- 
ated by §2. Since the curves for B are both above the lower 
bound, it is possible for the system to experience a decrease 
in fidelity under either type of jump. 



near 1/4. On the other hand, for other initial conditions, 
Ri for the /ii = U-^ solution is below that for the /ii = 1/2 
solution for all e < 1/4, and the same ordering as seen 
in the asymptotic limit (Fig[7|) occurs for some range of 
s near 1/4. Because — ln(C) diverges (slowly) to +oo as 
£ ^ for /ii = for a fixed initial condition this last 
solution always has the largest initial rate of convergence 
Ri diS e ^ 0. 

The last question we address in this section is the effect 
jumps have on fidelity. For the model we are considering 
we can easily evaluate the upper and lower bounds in 
Eq. (|48|) to be 



Amin = l/iiP + 1/2 - v/|MiP + 1/4 (59) 

Amax = l/ilP + 1/2 + + (60) 



For /ii = 1/2, these give the bounds 



0.043 : 



3 
4 




1.46. (61) 



To calculate B from Eq. ( [47| ), we observe that for /ii = 
1/2 we have observe that \0i\ = IO2I, \P2\ = |Mi| x 
This gives B = |/32p/|/3ip = 1/4. Because this is greater 
than 0.043 it follows that sometimes jumps can decrease 
the fidelity. For jui = f-^ we determine B and the bounds 
numerically. We find similar results, as plotted in Fig. |8] 
This time however B depends on which jump operator 
one is considering, and a fidelity-decreasing jump is much 
less likely following a stage of unstable evolution. Our 
results confirms what was seen in the second jump in 
the untypical evolution depicted in Fig. [3) where a jump 
following a stable stage of evolution visibly decreased the 
fidelity. 



C. Three-state jumping 

We now consider three state PR ensembles with cyclic 
jumps as in Ref. [61, still for a qubit subject to resonance 
fluorescence. In this case, in Eq. ([6| only Ki2^ /^23 and /^si 
are nonzero and Eqs. ([5|-(|6| yield a total of 12 equations 
(9 equations for jumping conditions and 3 equations for 
normalization condition) and 12 unknowns. Since this 
system involves quadratic equations, simple analytic so- 
lutions no longer exist. We find all three-state cycles by 
numerical search for all real solutions to Eqs. ([5|-(|6| 
using symbolic-numerical algorithms based on comput- 
ing a Groebner basis [16l [17], as described in detail in 
Appendix [BJ 

One of the reasons to study i^T-state PR ensembles for 
K > 2 \^ the search for low entropy solutions that allow 
for efficient tracking. As we have shown, for £ > 1/4 
there are no low entropy solutions ior K = 2. However 
by moving to = 3 we open the possibility for ensem- 
bles with different properties. The intuitive reason for 
the greater flexibility is as follows. Recall that for two- 
state jumping, we could use only one eigenvector of A 
to generate the PR states and we could not use complex 
eigenvectors because Bloch vectors must have real com- 
ponents. Observe that complex eigenvectors of A come 
in conjugate pairs. Thus, for three-state jumping, we can 
construct (real) Bloch vectors for PR states if we use con- 
jugate pairs of eigenvectors of A. This method imposes 
additional constraints however and cannot be used for all 
possible e. We now explore when complex eigenvectors 
can yield a PR ensemble. 

We first show that any cyclic jumps between three 
states is generated from two eigenvectors of A. We ob- 
serve that Eqs. ([6| imply that for Sj = Vj — rgg, with 
{1,2,3}, 

{A - K.i2){A - K.2'd){A - K3l)Sj = -A^12^23^3l5j (62) 

Since A is invertible, we can reformulate this condition as 
g{A)sj = 0, where we have defined a quadratic function 
g{A) = A^ - {ni2 + ^^23 + ^31)^ + {i^i2i^2?> + ^12^13 + 
^23^31 )• This is an eigenvalue equation for Sj. Note that 
all eigenvectors of A are also eigenvectors of g{A) and vice 
verse. Thus Sj is a linear combination of eigenvectors of 
A, whose eigenvalues A satisfy equation g{X) = 0. If all 
eigenvalues of A are distinct, then only two eigenvalues of 
A can satisfy equation ^(A) = and only two eigenvectors 
of A are used to construct Sj. 

However, not every conjugate pair of complex eigen- 
vectors of A can be used to construct Sj. Suppose A has 
two complex eigenvalues A and A*. Equation ^(A) = 
implies that /^i2 + /^23 + /^3i = 2ReA and 1^12 1^23 -\- 1^12 1^13 -\- 
^23^31 = |Ap. By construction, /^i2, 1^23^ ^31 are real, 
which is only possible if (ReA)^ > 3(ImA)^. In particular, 
for the resonance fluorescence example, complex eigen- 
vectors cannot be used to generate three-state PR states 
if \e\ > 1/2. This bound is actually not tight as was 
shown by numerical search for solutions to Eqs. ([5|-(|6| 
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FIG. 9: (Color online.) Solid arrows show Bloch vectors for 
3-state jumping. The volume of the sphere at the tip of each 
arrow represents the probability that the qubit occupies the 
corresponding pure state. The dashed arrow is fss- Solutions 
in (a) and (b) are shown for s = 0.27. These are low entropy 
solutions for 3-state jumping that are generated from complex 
conjugate pair of eigenvectors of A. Note that low entropy 
solution for 2-state jumping does not exist for such e. The 
solution in (a) has lower entropy than in (b). 



for three-state jumping [6j. These determined that such 
solutions exist iff < 0.282. 

Solutions for cyclic three-state jumping come in pairs. 
As \s\ approaches 0, the Shannon entropy h for half of 
the solutions approaches 1.2, whereas h for the other 
half approaches 0. For 0.247 < \s\ < 0.282, there are 
two preferred ensembles generated from complex eigen- 
vectors u± of A and shown in Fig. [9] with e = 0.27. In 
the region 0.183 < \e\ < 0.247, there are 6 solutions, 
which are shown in Fig. 10 with e = 0.23. Solutions 



with the same entropy are shown in one subplot. These 
entropy-degenerate solutions are mirror images with re- 
spect to X = plane and are constructed from ui and U- . 
The other two solutions have lowest and highest entropy 
and are generated by U-^ and U-. In the last region, 
l^l < 0.183, there are 8 solutions, which are shown in 



Fig. 11 with e = 0.18. Solutions with the same entropy 
still appear in the same subplot and are mirror images 
with respect to x = plane. All other solutions are con- 
structed from U-^ and U-. Solutions that have unique 
entropy always lie in x = plane. 

In ensemble {p/e, rk}k^i,2,3^ we can order the states 
by probability pk in decreasing order. The system al- 
ways jumps from one state to another in this order. For 
solutions with lower entropy, the systems spends most of 
the time in state one, which is nearly aligned with the 
steady state. For the lowest entropy ensemble, the states 
are spread out on the Bloch sphere as far away from each 
other as possible. As the entropy increases, states with 
small probability tend to move closer to the steady state. 
And for high entropy ensembles, the probability for occu- 
pying each state tends to equalize, and none of the states 
align with the steady state, but instead cluster around 
it. In Table [ij for £ = 0.15 we report the geometric char- 
acteristics of the six distinct solutions. The total angle 
ZSa = /(ri, + ^(^2, ^3) + ^(^3, n.) between the Bloch 
vectors show inverse correlation with the entropy. 





(c) 

FIG. 10: (Color online.) As in Fig. |9] but showing aU 3- 
state jumping solutions for e = 0.23. The solution in (a) has 
the smallest entropy. Every consequent solution has a larger 
entropy. Solutions (a) and (d) correspond to (a) and (b) in 
Fig. [9] Solutions (b) and (c) are actually pairs of solutions, 





/(ri,rss) 






h 


235.489 


115.323 


2.42085 


0.0313546 


0.020 


221.528 


109.238 


3.5805 


0.0390626 


0.023 


189.578 


94.7316 


5.87493 


0.0575965 


0.026 


26.9345 


0.654795 


12.8124 


5.30314 


0.466 


14.2435 


1.6635 


5.45759 


2.77901 


1.171 


13.0614 


3.06505 


1.4863 


3.46569 


1.299 



TABLE I: This table reports geometric characteristics of six 
distinct solutions for £ = 0.15. It shows that the total angle 
(column 1) for the ensemble shows the inverse correlation with 
the entropy h. 



D. Stability of three-state jumping 

We now analyze the stability of the above 3-state jump- 
ing schemes. We do this using results from Sees. |III C 



[III D| and numeric solutions for three-state jumping com- 
puted in Sec. |IV C[ As shown before, the stability of an 
individual trajectory is determined by the p rop erties of 
the effective Hamiltonian H^jii) from Eq. ([12|, where 
index i ranges from 1 to 3 and jii are the settings for 
adaptive monitoring that generate three- state jumping 
for resonance fluorescence system. We extract values for 
lii using the expression derived in Ref. [6] for cyclic jumps 
that states that 



(63) 



where /c+l stands for /c+l mod 3 and hk is some constant. 



Once \v^) and 



are known, coefficients ji^ and hk are 



uniquely determined. Thus our approach for determining 
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(b) 





0.00 0.05 0.10 0.15 0.20 0.25 

s 

FIG. 12: (Color online.) Ensemble Shannon entropy h for 
the eight different three- state jumping solutions (two of which 
come in degenerate pairs). The solid line indicate mean- 
square stability as well as piecewise deterministic stability. 
The dashed line indicates that the solution has only mean- 
square stability, with one of its stages being unstable. 





(f) 



FIG. 11: (Color online.) As in Fig. [To] but showing ah 3- 
state jumping solutions for s = 0.18. Again, the solution in 
(a) has the smallest entropy and entropy increases with every 
consequent solution. Solutions (a) and (b) correspond to (a) 
and (b) in Fig. 10 Solutions (e) and (f) correspond to (c) 
and (d) in Fig, 
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the stability of three-state jumping proceeds as follows. 
We convert the three- state jumping solutions in terms of 
Bloch vectors fi r2 and fs into state vectors 1'^^), l'^;!), 
and I'^l). Using Eq. (63), we compute /j^k and H{fik) 
with /c = 1, . . . , 3. At this point we can conclude that 
evolution stage from state \vl) to state is stable 

if state \vl) has the smallest eigenvalue with respect to 
operator H^jik)' Otherwise, such stage is unstable. From 
the eigenvalues A| and one can determine mean-square 
stability using Eq. (41). 



Results for this analysis are reported in Fig. [12] and 
resemble the results for two-state jumping quite closely. 



Figure 12 shows ensemble Shannon entropy h for all pos- 
sible three- state jumping solutions. Instability of indi- 
vidual trajectories for three-state jumping always arise 
in the same way as for two- state jumping: one of the 
three stages (corresponding to jumping from least prob- 
able state) becomes unstable. Solutions with small en- 
tropy (i.e., solutions with ^ as e ^ 0) have mean- 



square stability, but not piecewise deterministic stability. 
Solutions with large entropy have both mean-square and 
piecewise deterministic stability, except in small region, 
where the entropy curves approach lower entropy solu- 
tions and develop one unstable stage in the evolution. 



V. CONCLUSION 

In this paper we have considered a qubit undergoing 
evolution given by a Markovian master equation, subject 
to continuous monitoring that resolves every jump and 
allows the system to stay in a pure state. We studied 
special adaptive monitoring schemes that, in principle, 
allow an experimenter track the evolution of such system 
with finite-state machine as an apparatus. That is, the 
system jumps between only finitely many different states, 
the states in the associated physically realizable ensem- 
ble. The main contribution of this paper beyond Ref. [6] 
was analysis of the stability of such monitoring schemes. 
This is necessary to establish that the finite physically 
realizable ensembles introduced in Ref. [6 are not just 
mathematical constructions, but really are physically re- 
alizable. 

Because the evolution of the system is stochastic — it 
undergoes jumps at random times — we concentrated on 
the average properties of the system over many different 
realizations. Specifically, we derived conditions for the 
mean-square stability (that is, the long-time convergence 
of the fidelity of the system with the correct ensemble 
state). We showed that there is a positive parameter C, 
which is never more than unity, that guarantees mean- 
square stability as long it is not equal to unity. For the 
specific example of resonance fluorescence we considered 
11 different finite-state monitorings (with two or three 
different states). In all cases C was strictly less than 
unity (by a long way, in fact). Based on this, we con- 
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jecture that all finite physically realizable monitorings 
are mean-square stable. However we also considered the 
long-time rate R of convergence (which is exponential in 
time) and found that some monitoring schemes converge 
much faster than others, in a way contrary to what is 
suggested by considering C alone. 

Although all monitorings we considered were mean- 
square stable, there is a variety of behavior in the conver- 
gence of individual trajectories. Some monitorings give 
trajectories that are guaranteed to be piecewise deter- 
ministically stable (i.e. the between-jump evolution al- 
ways increases the fidelity). However others do not — 
in some deterministic (between jump) stages of the evo- 
lution the system moves away from the ensemble state 
that it "should" be in at that stage (and will be in, in 
the long time limit). Moreover, the stochastic jumps can 
contribute to both stability and instability of individual 
trajectories, even under piecewise deterministically sta- 
ble monitorings. However in the long time limit, as the 
system approaches the ideal ensemble states in the mean- 
square sense, a fidelity decreases upon a jump become 
very unlikely. 

We illustrated these effects for the specific example 
of resonance fluorescence with an adaptively controlled 
weak local oscillator. Interestingly we found that those 
monitorings that are not piecewise deterministically sta- 
ble tend to be those that are most stable in the mean- 
square sense (by the measure of C being very small). 
Moreover, these are the monitorings that produce low en- 
tropy solutions, where the system spends most of its time 
in one state in the ensemble. Such monitorings would 
make it possible, in principle, to track a large number N 
of qubits, using much less than N bits of memory (i.e. a 
finite state machine with far fewer than 2^ states). 

It is important to note that there are many open ques- 
tions in this field of quantum state tracking with finite 
state machines. First, given a physically realizable finite 
ensemble, does there exist an explicit construction for an 
adaptive monitoring scheme that realizes it? Second, is 
any such monitoring mean-square stable, as conjectured 
here? Third, is it true that any D-dimensional ergodic 
Markovian quantum system can be tracked by a K-state 
machine with some finite i^>(I) — 1)^ + 1, as conjec- 
tured in Ref. [6l. Fourth, is = (D - 1)^ + 1 always 
sufficient? Fifth, can an example system be found that 
proves that K = D is not sufficient in general {D = 3 
would be the minimum system size for such a search). If 
the last can be proven, then there may be a relation to 
the recent result that there are classical stochastic pro- 
cesses that can be generated using quantum systems of 
lower entropy than that required using only classical sys- 
tems [T9^. Finally, it seems likely that any given master 
equation would have additional structure that would en- 
able one to use a smaller K than that conjectured above, 
and this idea also remains to be explored. 
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Appendix A: Useful Identities 

We now explain how to derive expressions for the 
convergence coefficient C that determines how quickly 
the state subject to adaptive monitoring approaches the 
states in the corresponding PR ensemble. We begin with 
the two-state jumping scenario. In this case, the expres- 
sion for C is given in Eq. p3|. We evaluate it through 
relating QliQii and Q22Q22' We do this by comparing 
the general expression J d^r\/3n\'^ to J d^r\an\'^ for the 
special case when the initial state |?/^o) is the jumping 



stat \vi)^ i.e. = 1 and /3o = 0. Then with Eqs. (23) 



(25) and Eq. (15), we get 



0^21+1 = Oi2ie 



-Air2H 



Oi2l = Oi2l-ie 



-^2^21 



Again by induction we deduce that 

I 

a2i = iQl,Ql,yi[e-^'^-^^-^e-^'^-^^ 



a2i^i=Qli{QliQliyile- 



(Al) 
(A2) 



(A3) 
(A4) 



Then 



By construction j d?^T\\\ilj2i) 
proves that 

\Q\iQ\i\ 



|2 = 1. Therefore, Eq. (AS) 



4ReAfReA| 



1. 



(A6) 



To calculate the value of C, we relate QiiQii to (322^22- 
We recall the discussion after Eq. (17), which showed that 



operator 52^1 is proportional to the identity operator. 
Using this fact and Eqs. (25), we see that 



S2SM) = QliQliK) 

S2S1K) = gLQ22l<)- 



(A7) 
(A8) 



Since 52^1 is proportional to the identity, we can conclude 



?iiQii 



^22^22- 



This discussion can be easily generalized to the in- 
state jumping scenario. We again calculate the value 
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of / d^^T\anK\'^ in two ways when ao = 1 and /3o = 
0. The first way is to note that, as before, in this case 
/(i^^r|anKp = J d''^r\\\ijjnK)\\^ = 1- The second way 
is to use the fact that in this case 



n-l 



anK = {Qli ■ ■ ■ gfi)" n • • • e-^-^< 

so that 



/ 



jnK-?i |2 _ / Qll * ' ' Qn 



2^ReXl ■ • • ReA|^ 



(A9) 



(AlO) 



Thus Qn---Qfi = 2^ReAf •••ReA|^. Next, we show 
that Qli • • • = Q22 ' • • Q22' We prove this relation- 
ship for Ql.j^ by considering full cycle jump operator 

Sk = sk - • -si. 

Assuming (as we can do without loss of generality) that 
the jump operator c is traceless, we have cx / so that 

SkK) = [^i(/ii,...,/iK)+^2((/ii,...,/iK)c]|vi) (All) 

The exact nature of functions gi and ^2 depend on the 
proportionality constant relating to the identity oper- 
ator. Just as for Eq. ( pT] ), a system can undergo jumping 
dynamics iff ^2 = 0- This means that the full cycle jump 
operator Sk is proportional to the identity and we con- 
clude that Oil • • • Qn = Q22 ' ' ' using the trick from 



Eqs. (A7)-(A8). 



Appendix B: Groebner Basis 

The task of finding jumping states for adaptive un- 
ravellings requires solving non-linear equations. We now 
introduce important concepts from computational alge- 
braic geometry and review some algorithms used to solve 
multivariate polynomial systems. Our presentation will 
rely on analogies with linear algebra. From now on we 
regard a polynomial as the finite sum of terms, where 
each term is a product of a coefficient and monomial. 

Suppose we want to solve the system of non-linear 
equations, 

{h{Xn) = 0, f2{Xn) = 0, . . . , /,(x^) = 0}, (Bl) 

where {/i, /2, . . . , /s} are the polynomials with real ra- 
tional coefficients and Xn = (xi, X2, . . . , x^) is the list of 
variables. Our goal is to find all the solution to the set 
of equations in Eq. ( |B1[ ) and a concept of ideal becomes 
useful. A collection of polynomials generates the ideal via 

I = (/l, /2, . . . , /n) = {ELl ^ifi ' hi,...,hs e C[fn]}, 

where s can be any finite index and C[xn] is a collection 
of all possible polynomials with complex coefficients with 
variables Xn- Thus ideals are similar to vector spaces, 
which are formed from all possible scalar combination 
of vectors. But instead of scalars, one uses all possible 
polynomial functions, hk^ defined on C[xn] to form an 



ideal. Ideals are important because the solution set to 
the newly created ideal and to the original system are 
the same. 

Different set of equations can have the same solution 
set and thus generate the same ideal. One of the main 
ideas in algebraic geometry is to pick a good set repre- 
senting the ideal that has nice properties and yields easy 
way for identifying the solutions. Such a set is called a 
Groebner basis. 

Gaussian elimination is the algorithm used to solve sys- 
tem of linear equations. The extension to this algorithm 
used to solve a system of polynomial equations is known 
as Buchberger's algorithm and it is implemented in many 
packages for symbolic computations such as Mathemat- 
ica. Maple, Sage. The set of equations obtained as a 
result of these algorithms is known as the Groebner ba- 
sis. 

Before proceeding we explain how to check if a given set 
of equations produced by some software package is indeed 
a Groebner basis. To do this, we first need to know the 
S-polynomial. Given two polynomials / and let be 
the least common multiple of leading terms of / and g, 
denoted by LT{f) and LT{g). Then the S'-polynomial 
is computed via S{f,g) = f^'^/LT{f) - g^'^/LT{g). A 
set of polynomials G = {^i, . . . , are a Groebner basis 
iff for all i ^ j, the remainder on division of S{gi,gj) 
by G is zero. An alternative explanation proceeds as 
follows. In general, one can always write S{gi^gj) as 
^{9ii9j) = Xir=i ^i9i~^^^ where are some polynomials 
and degree of the remainder polynomial r is smaller than 
any polynomial in G. Then G is a Groebner basis iff the 
remainder polynomial r is zero for all i j. Division by 
a collection of polynomials is usually implemented in the 
software packages for computational algebraic geometry. 

Now we explain how to use a Groebner basis to find 
all solutions to the system of non-linear equations. Un- 
like linear algebra, where a vector space always has the 
same number of basis vectors, different Groebner bases 
can have different number of elements and drastically 
different properties for the same ideal. The collection of 
Groebner bases arise from different orderings of monomi- 
als in the system of polynomials. Two particular order- 
ings are relevant for the current discussion: lexicographic 
order (Lex) and degree reverse lexicographic order (DRL). 
Lexicographic order is an alphabetical order (write out 
monomial in full without any powers and order like the 
words in a dictionary from left to right). For DRL order, 
we first compare the total degree and then perform lexi- 
cographic ordering by reading the expressions from right 
to left. Both orderings allow different access to informa- 
tion about solutions to polynomial system of equations. 
The details on different orderings can be found in [16] 

Lex order is particularly useful because of the elimi- 
nation theorem [T6j, which states that the Groebner ba- 
sis computed with respect to Lex order will yield a set 
of polynomials that can be solved by back substitution. 
This means that if the system of non-linear equations 
has finite number of solutions and the Lex Groebner ba- 
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sis are given by G = then gi = gi{xi) is 

a polynomial function of one variable only and we can 
determine all values for xi by solving gi{xi) = 0. The 
next polynomial is most likely ^2(^1,^2) is a function of 
two variables so that we can solve for X2 given all possible 
values for xi. This sequential substitution allows to solve 
for all variables x^- Some of the elements of G are just 
constraints that eliminate possible solutions. Such pro- 
cedure is used to solve a system of non-linear equations 
in Mathematica (command NSolve). 

This algorithm comes with some complications. For 
some systems, algorithm becomes unstable. This is so 
because the degree and size of coefficients in the Lex 
Groebner basis quickly become huge even for relatively 
small systems. And solving gi{xi) = involves rounding, 
which introduces errors that quickly propagate through 
back substitution and yield wrong results. 

The Groebner basis computed with respect to DRL 
ordering offers an alternative way to compute solutions 
to polynomial system of equations that does not involve 
back substitution. We now review machinery from alge- 
braic geometry needed to compute solutions using DRL 
ordering of monomials. The ffist construct involves the 
space of all polynomials C[xn] and ideal / defined above. 
Polynomials C[xn] can be classified into distinct cate- 
gories (cosets) such that polynomials / and g belong to 
same coset iff /— ^ have remainder with respect to ideal 
/ or equivalently polynomials / and g have the same re- 
mainder with respect to /. We denote the reminder of / 
by [/] . Note that for checking this criteria it is sufficient 
just to consider remainder with respect to the Groebner 
basis (same ordering as the one used for division algo- 
rithm) spanning /. 

We can consider the space of cosets, i.e., the space of 
the remainders with respect to division by /. This is a 
quotient ring, denoted by ^ = C [xn]- It also happens to 
be an algebra that has vector space structure and carries 
important information about solutions to the system of 
polynomials that spanned the ideal /. 

We now point the key steps in computing solutions to 
polynomials {/i, /2, • • • , /n} and illustrate them using a 
simple example. 

• We compute the Groebner basis, G = {^1, . . . 
with respect to DRL ordering for the ideal / = 

(/l,/2, • • • ,/n)- 

• We identify leading term for every element of the 
Groebner basis LT{gi). 

• We create set B from monomials that a not divis- 
ible by LT{gi)^ i.e., monomials in B have degree 
smaller than the degree of LT{gi). This set is a 
basis for space of remainders, A. 

For example, suppose we want to solve a system of equa- 
tions 

h=x^ -y^ ^xy = {) 
/2=x27/ + ^-l = 



The Groebner basis (DRL ordering and x > y) for this 
system is 

gi=x'^ ^xy- 

g2=y^ -xy'^ + y-l (B3) 
^3 = / + + 27/2 - X - 2y 

Leading terms for these polynomials (with respect to 
chosen ordering) are LT{I) = {x'^ ^xy'^ ^y^} and the set 
B = {1, X, 7/, XT/, 7/2, 7/^} is the basis for Quotient ring. 

To extract the information about solutions from the 
quotient ring we associate each polynomial / in C[xn] 
with a linear map rrif : A ^ A whose action is given by 
mf{g) = [fg]^ where g is some polynomial. Map m/ can 
be represented as a matrix. To do this, we consider the 
action of m/ on basis elements in the set B. And mf{h) 
for every & G B is an element of quotient ring A can be 
written as a column vector with respect to B. 

For example discussed above we show how to compute 
nix. The action of rrix on 1 is given by 

m,(l) = [x] = (0,1,0,0,0,0)^, (B4) 

i.e., rrix is the coset [x], which can be represented as a 
vector with respect to elements in basis set B. Map 
will take some basis elements outside of B. In this case, 
we compute the remainder with respect to the Groebner 
basis and again compute a column vector with respect to 
B. For example, 

mx{x) = [^2] = [g^^y'^-xy] = [y^-xy] = (0,0,0,-1,1,0)^. 

(B5) 

Then the matrix form for is 



(0 
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-1 


o\ 


1 














1 











-1 
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1 





-1 


1 








-1 





1 











-1 


Vo 








1 


1 


0^ 



(B6) 



We need to compute these matrices because they have 
very nice property: eigenvalues of m/ are the values of 
polynomials / on solutions for original solutions. Thus 
eigenvalues of rux are all possible values a variable x takes 
in the solutions set. Moreover, we don't need to compute 
every possible matrix m/ because rufg = mfrrig. This is 
significant since computing ruf relies on costly symbolic 
calculations, which are slow, whereas matrix multiplica- 
tion is fast. 

Thus eigenvalues for all different maps . tells us all 
possible values each variable Xi can take. But to learn the 
full solution (how to combine different Xi to form Xn) we 
also need information that is encoded in the eigenvector 
of m/ for some /. Let v be such an eigenvector and we 
normalize this vector so that the first component is one 
(assuming that 1 is the first element in B). If variable Xi 
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is the j-th element in the set then jth component of v 
is the value of Xi that yields one of the solutions. All val- 
ues extracted in this way from one eigenvector come from 
one solution. However, this procedure does not work for 
every polynomial /. If m/ has degenerate eigenvalues 
then corresponding eigenvectors cannot be used to ex- 
tract the solution. However, generic linear combination 
of the variables will yield a desired matrix. The details 
on this can be found in [171 HI] 

For three state jumping, we computed the solutions 
using the Groebner basis with DRL order. The set of 
equations for cyclic three- state jumping has a lot of sym- 
metry, e.g., one can map r i r 2 ^ r 3 ^ fi and 



K,i2 ^ K,23 ^ f^is ^ f^i2 and the system will remain un- 
changed. As the result, all matrices m^^. associated with 
unknowns in the system have degeneracy with respect to 
real eigenvectors and cannot be used to extract the solu- 
tion. Instead we used matrix m/, where / = rn + hzi2. 
Here rn is the first component of Bloch vector fi. The 
reason for this choice is the following. Matrices asso- 
ciated with f<ii2 and ^^23 and ^^31 are nice because they 
have least non-zero elements which speeds up the calcu- 
lation (symbolic part). However, any combination of K:12 
and i<i23 and ^^31 or same coordinate of ri r2, does 
not break the degeneracy due to symmetry above, but 
ru + ^12 does. 
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